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Abstract 

This paper is one of a series of papers describing the development of a new numerical 
approach for solving the steady Navier-Stokes equations. The key features in the current 
development are (i) the discrete representation of the dependent variables by way of high 
order polynomial expansions, (ii) the retention of all derivatives in the expansions as un- 
knowns to be explicitly solved for, (iii) the automatic balancing of fluxes at cell interfaces, 
and (iv) the discrete simulation of both the integral and differential forms of the governing 
equations. The main purpose of this paper is, first, to provide a systematic and rigorous 
derivation of the conditions that are used to simulate the differential form of the Navier- 
Stokes equations, and second, to extend our previously-presented internal flow scheme to 
external flows and nonuniform grids. Numerical results are presented for high-Reynolds- 
number flow (Re = 100,000) around a finite flat plate, and detailed comparisons are made 
with the Blasius flat plate solution and Goldstein wake solution. It is shown that the error 
in the streamwise velocity decreases like r a Ay I. 2 , where a ~ 0.25 and r = is the grid 
aspect ratio. 


I. Introduction 

This paper is concerned with the continued development of a new numerical approach 
for solving the steady Navier-Stokes equations and other conservation laws of the form 


dh x dh* 
dx "I" dy 


= 0 . 


( 1 . 1 ) 


The key features in the current development are (i) the representation of the discrete 
unknowns by way of high order polynomial expansions; (ii) the retention of all derivatives 
in the polynomial expansions as unknowns to be explicitly solved for; (iii) the automatic 
balancing of fluxes at cell interfaces without the use of flux limiters or reconstruction; and 
(iv) the discrete simulation of both the integral and differential forms of the governing 
conservation laws. 

In a previous paper (Scott et al, 1995), we formulated a quadratic-expansion discretiza- 
tion for the two-dimensional, compressible Navier-Stokes equations. A specific scheme was 
constructed for laminar channel flow, and it was shown that the developing boundary 
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layer could be accurately resolved on relatively coarse grids (as few as six cells per channel 
width). 

The main purpose of this paper is, first, to provide a systematic and rigorous deriva- 
tion of the conditions that axe used to simulate the differential form of the Navier-Stokes 
equations, and second, to extend our previously-presented internal flow scheme to external 
flows and nonuniform grids. Specifically, we (i) derive error bounds for the order- N Taylor 
series expansion of the exact solution to equation (1.1); (ii) formalize what is meant by 
a discrete polynomial solution to a conservation law and introduce a generalized concept 
of convergence; (iii) derive necessary conditions for a discrete polynomial solution of (1.1) 
to converge to the corresponding exact solution; and (iv) construct a scheme for high- 
Reynolds-number flow past a flat plate airfoil, and show through comparison with the 
Blasius solution that the streamwise velocity is obtained to order r“ Ay 2 , where a « 0.25 
and r = ^ is the grid aspect ratio. 

Sections II A. and B. below serve primarily as a review of our original formulation 
(Scott et al, 1995). However, here we limit ourselves to incompressible flow, and introduce 
a generalized flux vector notation which simplifies some of the derivations. In Section II 
C., we deal with items (i) - (iii) as outlined in the preceding paragraph, and Section III 
deals with item (iv). 


II. Discretization of the Navier-Stokes Equations 


A. Incompressible Navier-Stokes Equations 

We consider the two-dimensional, steady, incompressible Navier-Stokes equations in 
dimensionless form. We assume that the viscosity p is constant, and denote the density 
by p and the Reynolds number by Rcl, where Rzl — pu ™ L . The parameters L and 
refer to some reference length and velocity, respectively. 

Let x and y denote the horizontal and vertical coordinates, respectively, of a two- 
dimensional Euclidean space £ 2 - Denoting the horizontal velocity component by u, the 
vertical velocity component by v, and the static pressure by p, the governing equations 
for the conservation of mass and momentum may be written in Cartesian coordinates as 
(Anderson et al 1984) 

fht rh ) 

= 0 


du 

dx 


+ 


dv 

dy 


d ( 2 

al (u 


+ P - r xx ) + 


— (uv - T XV ) = 0 


d , 

dl (uv - 


a 


'xj 


where 


) + + P ~ Tyy) = 0 


2 ( n du _ dv, 
3Re L ' dx dy 


( 2 . 1 ) 


( 2 . 2 ) 


(2.3) 


(2.4) 
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(2.5) 


— 1 ,0u 

Tty = ikl c a^ + 0^ 

° Y 2 — - — ) 
L { dy dx ] 


Tyy ~ 3 Re 


( 2 . 6 ) 


(Although || or may be eliminated from t xx and r yy using (2.1), we retain both terms 
here to be consistent with our original compressible formulation.) 

The integral form of equations (2.1) - (2.3) is given by 


l 

l 

i 


S(V) 


S(V) 


o 

II 

t-3 

* 

(2.7) 

—* ■■ ■■ fr 

h X M 'ds = 0 

(2.8) 

h y M * ^ ~~~~ 0 

(2.9) 


S(V) 

where the flux vectors, h M , h XM , and h YM , corresponding to the conservation of mass, 
x-momentum, and y-momentum, respectively, are defined by 


h M = f (u,v) 


( 2 . 10 ) 


h X M d — (« 2 + P - T xx ,uv - T xy ) (2-H) 

h YM d = (utl — T xy , V 2 + P — Tyy). (2-12) 

S(V) is the boundary of an arbitrary region V in E%, and ds is equal to da n, where n 

is the outward uni t normal to S(V) and da is the length of a surface element of S( V ) . 

Equations (2.1) - (2.3) and (2.7) - (2.9) are the differential and integral conservation laws, 
respectively, for the conservation of mass and momentum in two space dimensions. 


B. Discrete Flux Conservation Equations - Integral Formulation 

Let E 2 be discretized by a mesh with nonoverlapping rectangular regions. We assume 
for the time being that the mesh spacing is uniform in each coordinate direction. (See 
Figure 1.) Each of the rectangular regions in the mesh is referred to as both a conservation 
element and a solution element (Chang, 1995). A conservation element is a discrete region 
in E 2 over which the discrete analogue of the integral conservation laws (2.7) - (2.9) 
is imposed. A solution element is a discrete region in E 2 in which a local polynomial 
expansion is employed to represent the physical solution. In general, they need not refer 
to the same discrete region. A conservation element is denoted by CE(i,j) and a solution 
element by SE(i,j). The boundary of a conservation element is denoted by S(CE(i,j)), 
and its cell center by (x;,yj). 
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We assume that the velocity and pressure can be represented on each SE(i,j) by a 
polynomial expansion of degree two about the cell center (x,-, yj), as follows: 

«(*>y;*\j) d = uo.o + + «o,i(y-yi) (2.13) 

+ u 2 ' 0 (x-Xi ) 2 + u ltl (x — Xi)(y — yj) + u 0 i 3 (y-yj ) 2 

y(x,y-,i,j ) = v 0 ,o + r 1>0 (x-x,) + v 0il (y-yj) (2.14) 

+ Vs,o(x - Xif + v Xil (x - x,)(y - yj) + t> 0 , s (y - yj) 2 


p(x,y;i,j ) = Po.o + Pi t o(x-xj) + p 0 ,i(y-y>) 


(2.15) 


+ Ps,o(* - Xi) 2 + p 1>x (x - x<)(y - yj) + p 0 ,s(y - yj ) 2 . 


For clarity, the t,j subscripts have been omitted from the coefficients in the local poly- 
nomial expansions. These coefficients Eire the unknowns to be solved for. They are related 
to the discrete derivatives (at the cell center) by 


u 

(2.16) 

du/dx 

(2.17) 

du/dy 

(2.18) 

^ d?u/dx 2 

(2.19) 

d^u/dxdy 

(2.20) 

\ d?u/dy 2 , 

(2.21) 


and similarly for v and p. 

The discrete analogue to equations (2.7) - (2.9) in E? is then given by 


<p h M • ds = 0 

Js(CE(i,j)) ~ 

(2.22) 

f 

y h X M • ds = 0 

S(CE(i,j)) " 

(2.23) 

r ^ ► 

p h Y m • ds = 0 

S(C£(»,i)) ' 

(2.24) 
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where 


h M d = («,v) 


hx m — (*j "I - P T* x » If ~ T*y) 


7* dtf / 2 . \ 

h Y M = (*}tj r*y> y + p ) 


and 


7 * -r <r — 


3 Rgl 

1 


1151 Re t 


(2 du/dx — dv/dy) 
-( du/dy + dv/dx ) 


Ivy = 


ZRe L 


(2 dv/dy - du/dx). 


(2.25) 

(2.26) 

(2.27) 

(2.28) 

(2.29) 

(2.30) 


Equations (2.22) - (2.24) are a coupled system of integral conservation laws in which 
the fluxes h M • ds, h XM • ds, and h YM ■ ds are conserved by way of the discrete variables 
u, v and p. Each equation takes the form 


i 


S(CE(i,j)) 


h ■ ds = 0, 


(2.31) 


where the second-order expansion h is a function of u, ij and p. Since the form of 
each integral in (2.22) - (2.24) is identical, the integrations may be carried out by way of 
equation (2.31), where 

(2.32) 


h = f (h x ,h y ) 


and 


h x (x,y-,i,j ) d = h x 0 + h* 0 (x-Xi) + K.iiy-Vj) (2.33) 

+ h x 0 (x — z,) 2 + h x tl (x - Xi)(y - yj) + h x 3 (y-yj) 2 


h y (x,y;i,j) = f h y i0 + h?. 0 (*-x;) + h'^y-yj) (2.34) 

+ h y 2 0 (x-xi) 2 + h y t l (x-xi)(y-yj) + h y 2 (y - yj ) 2 . 

It is understood that each of the coefficients in (2.33) and (2.34) arejunctions of the 
discrete variables u 0 , 0 , v 0 , 0 , p 0 , o, «i,o, v 1>0 , Pi.o, etc.. For example, when h corresponds to 
h XM , the term corresponds to the constant term of the expression u 2 + p — t xx , and 
similarly for the other terms. Once the results are obtained in terms of h, it is a simple 
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matter to derive the corresponding results for h M , h xti , and h YM by way of equations 
(2.25) - (2.27). 

Since the boundary S(CE(i,j )) of each conservation element is a simple closed curve 
in E 2 , the surface integration required in equation (2.31) can be converted into a simple 

contour integral. With ds — daft, where n is the outward unit normal to S(CE(i,j)) 
and da is the length of a surface element in E 2 , we have 


ds d = dy i — dx j , . 


(2.35) 

(2.36) 

(2.37) 

(2.38) 

The line integration is taken to be positive in the counterclockwise sense. If we denote 
the vertices of an arbitrary conservation element CE(i,j ) by P, Q, R, and S as shown in 
Figure 2, we have 

<p h • ds = <p g • dr 

Js(CE(i,j)) ' JPQRSij ~ 


and 


where 


and 


h- ds = — h y dx + h x dy = g ■ dr 


9 = (- h*,h*) 


dr d = dx i + dy j. 


d = f [j(PQ) + J(QR ) + J(RS) + J(SP)]_. 


(2.39) 


[J(PQ)]i,j denotes the flux of h through the line segment PQij, and similarly for J(QR), 
J(RS ), and J(SP). We then have (omitting i,j subscripts) 


rQ 


Similarly, 


J(QR ) = J 


h v dx 

r x <+^r 

+ h x dy = / 

J ij — ^ 

h y dx 

with y = 

def 

/•yi + T 1 

- h x dy 

with 

X X j 

def 

— / h 9 dx 

with 

1 

II 


J(SP ) d = 


ry>+~ 3 1 

/ V 


dy 


with x = Xi + 


Ax 

~ 2 ~ 

Ay 

2 

Ax 

~ 2 ~' 


(2.41) 

(2.42) 

(2.43) 
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Carrying out the line integrations in equations (2.40) - (2.43), one obtains 

J(PQ) = ^~ h U + Ax [^f- h’, + ^A>, + A» 0 ] 


J(QR) = - ^ K, - Ay K. - X + 


Ax 


Aar 3 

J(i!S) = - =r- hi 


12 ..... - Ax [^fhl, - ^AJ, + fc»„] 
J(3F) = + Ay ^ kl . + A?,.]. 


By virtue of equations (2.31) and (2.39) we require that 

J(PQ ) + J(QR ) + J(RS) + J(SP) = 0. 

Thus, we obtain the flux conservation constraint 

hU + hl A = 0. 


(2.44) 

(2.45) 

(2.46) 

(2.47) 

(2.48) 

(2.49) 


Imposing this condition, we obtain the following expressions for the normalized flux of h 
across the boundaries of CE(i,j ): 


J(PQ) __ Ax 2 Ay 2 Ay x 

“aT" _ TF + 4 2 + 


_ J(qg) 

Ay 

_ JQRS) 
Ax 

J(5P) 


12 

Ay 2 


4 
Ax 2 


12 + 4 '‘ s>0 


a;. - + K, 


A x 2 iy . Ay 2 , y Ay , x l y 

: 12 ^i.o + 4 ^ 0,2 + 2 ^ 1>0 

Ay 2 , r Ax 2 x Ax x x 

*^ 0,2 4 ~ ^ ** 2,0 4 “ q « 1,0 ■ " 0 , 0 - 


(2.50) 

(2.51) 

(2.52) 

(2.53) 


Ay T ‘ 12 "°* s ' 4 " so ‘ 2 

The flux expressions above may now be expressed in terms of the discrete dependent 
variables of the Navier-Stokes equations by way of equations (2.25) - (2.27). Corresponding 
to (2.25), we have 

- - (2-54) 


h — hu 


so that 


h x = u h v = v 

and the mass flux conservation constraint corresponding to (2.49) is 

tii.o + U 0 ,1 = 0. 


(2.55a, b) 


(2.56) 
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The normalized fluxes for h M are given by (2.50) - (2.53), where h* 0 = u 0i0 , ho.o = v 0 , 0 , 
etc. 

Corresponding to (2.26), we have 


h = h 


•iXM 


(2.57) 


where 


h x = u 2 + p — r a 


h* = uv — r zy . 


(2.58a, b) 


The x-momentum flux conservation constraint corresponding to (2.49) is 


1 2 

«i,o«o,o + «o.i u o t o + Pi,o - -^-[(2u 0 ,2 + t>i.i) + 3 (4u 2> „ - v l>x )] = 0, (2.59) 

and the flux expressions for h xu are readily obtained by way of equations (2.58) and 
(2.50) - (2.53). Thus, one easily obtains the x-momentum flux through the right hand face 
of a conservation element, _____ 

J(SP) XM 


Ay 

-^|~(2u 0iJ u 0> o + «0,1 + Po.i) + — (2u 2>0 u 0f0 + u* 0 + p 2i0 ) 


(2.60) 


- ^ [fo.o u ° (l + tio.oUo.i - — ^-( 2 u 0 , 2 + U 1(1 )] + U ».« + Po.® ~ 3^(2^..-^), 

and similar expressions for the remaining fluxes. (Here A* # has been replaced with — ho.i 
because it gives a simpler expression.) 

Similarly, corresponding to (2.27), 



h Y u 


(2.61) 


with 


h x = uv — T 


L*V 


h* = v 2 + P - In- 


(2.62a, b) 


The y-momentum flux conservation constraint is 

1 2 

t>i,o«o,o + v 0il v 0i0 + p 0il - -^-[(2 v 2i0 + Ui,i) + 3 (4v 0>2 - u ljX )] = 0, (2.63) 

and the flux of h YM through the right hand face of a conservation element is given by 


J(SP) 

Ay 


YM 


(2.64) 
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Ay 2 

12 


Ax 2 


(u 0 ,2«0.0 + «0,!«0,0 + Uo.lUo,!) ■+ 4— (t>a,o«o,o + Mj.oUo.o + Vi.ott lf o) 

+ [ U0, ° Ul ’° U °-° Ul ' 0 “ ^ _ ( 2Vl -° +Ul - 1 )] + «0, V „.o - +«!,«)• 

The following conditions are then satisfied on each conservation element: 


J(PQ)m + J{QR)m + J(RS)u + J(SP)m = 0 (2.65) 

J(PQ)xm + J(QR)xm + J(RS)xm + J{SP)xm = 0 (2.66) 

J(PQ) ym + J(QR)ym + J(RS) ym + J{SP)ym = 0- (2.67) 

The above formulation provides the framework through which integral flux conserva- 

tion is achieved. We may now turn our attention to the differential conservation laws (2.1) 

- (2.3). 


C. Discrete Flux Conservation Equations — Differential Formulation 
In the previous section, we derived the local constraint (2.49), 


K. + *■?.. = °. 

to ensure that the discrete flux vector h (defined by (2.32) - (2.34)) satisfies 


/ 

Js(CE(i,j ;)) 


h - ds = 0. 


Now we observe that, by virtue of (2.49), h also satisfies 


V • h = 0 at (x,y) = (xi,yj), 

i.e., the differential conservation law is satisfied at the cell center. 

Suppose that we further require h to satisfy V • h = 0 identically throughout each 
SE(iJ). Then one obtains two additional local constraints, 

2h* 0 + hi, = 0 (2.68) 

+ 2 h* 0>3 = 0. (2-69) 

Our main objective in this section is to rigorously derive the above two constraints 
from first principles. We show, in particular, that enforcement of these constraints ensures 
the satisfaction of necessary conditions for the discrete second derivatives to converge to 
their analytical counterparts. 
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Throughout this section, we shall make use of the following result, which may be 
established by induction. 

Lemma 2.1 Let h = (/i z , h 9 ) be defined and continuous on an open domain V of 
E 2 , and let the partial derivatives of h x and h 9 exist to all orders and be continuous on 
V. Then, h satisfies V -h = 0 throughout V if and only if its partial derivatives satisfy 


8 n h r d n h? 

ax .-»a yt (»-») + 3l n-,-, v+ i (^y) = 0 

for n = 1, 2, 3, ... and k = 0,1, 2, ..., n — 1, 


(2.70) 


for every (x, y) in T>. 

As a preliminary to what follows, we first consider the N - th order Taylor series ex- 
pansion of the exact solution h of the conservation law 


V-H = 0. 


(2.71) 


We assume that h — (h x ,h 9 ) satisfies (2.71) throughout some open domain V of E 2 , and 
that h is analytic throughout T>. Let h N denote the JV-th order Taylor series expansion of 
h about a point (x 0 ,y 0 ) in L>. Then 


lx,_ \ dnh * < \1 ( x ~ x ») n * (y “ Vo) 1 

M*>y) - (n — k)\ k\ 

n=0 Jfc=0 


and 


Since 


and 


UV( \ - r ( .1 (x-x„) n k (y — y o y 

K( x ,y) zJZ,[a x n-* 5y *( x ®’^)J ( n _ fc)i fci 

n= 0 k = 0 v 7 


dh JL - VV'f dn b X 


dx ^ 4 Ldx n ~ k dy k v ’ (n — A; — 1)! 

n= 1 k = 0 v 7 


( z (y-y.) 

1 ,,V "J (n-k-lY. Jt! 


dh y N 

dy 


N n 


^ f d n h* _ ,1 (x - x 0 ) — ty - y.)'- 

- ^^Ldx»- fc ay^ 0,y ° ; J (n-Jfc)! (k — 1)! 


(X - X.)"“* (y - y«) (t_1) 


N n— 1 


d n h y ..Q fx-xo y (y-y 0 r 

- 2^ Z^[dx n - k - l dv k + lK °' yo) \ (n — k — 1)! k\ ' 


(x - xo)"-*- 1 (y - y 0 )* 


n=l fc=0 


dy*+l V-o, »o,j 


(2.72) 


(2.73) 


(2.74) 


(2.75) 


it follows from Lemma 2.1 that h N is a solution of equation (2.71). We will refer to h N as 
the exact order-N solution of (2.71). 
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Let (x, y) be another point in V. Then, as an approximation to the exact solution h 
at (x,y), the error in h N is given by 

, , (x - x 0 )* +1 -* (y - yo) fc . . 

h x (x,y) - h N (x,y ) - (jv + 1 - k)l k\ 

d” +1 h* , . , xl (x - x 0 )" +1 -* (y - y«)* 

h*(x,y) - h* N (x,y) = EU x *+i-*^ (x 2’M (jv + i - jfc)T Jfc! ( ??) 

where (xj , yj) and (x^yj) are points on the line segment between (x„, y 0 ) and (x, y) (Buck, 

1978). 

Now suppose that 


d N+1 h x 


\dx N + 1 ~ k dy k 


< Mf,+i t k and 


d N+1 h y 


dx NJrl ~ k dy k 


< Mff+i'k 


(2.78) 


in a neighborhood Nq of (x„, y 0 ) for k — 0,1, ..., JV +1. If Ax — (x— x 0 ) and Ay — (y y 0 ) 
where (x,y) is any point in No, then 

li? ? || < VM |A*r*-* |Ay| fc 

||/i M°o < 2^ M ” + i> fc (jv + l-Jb)! Jk! ' 


(2.79) 


If we let 


M N+l — suplA/^+iji fc — 0, 1, ..., .AT 1}, 
we have the more conservative error estimate 

M r r Mn+i (N ± 2) [max([Ax|, |Ay|)] N+1 

II" "jNr||oo < 


(2.80) 


(2.81) 


where [(A±±)!] 2 d = (f )! (^)! if N is even. 

The exact order-JV solution h N is the prototype discrete polynomial approximation 
to the exact solution h. Its order of accuracy is given by (2.81). We may now turn our 
attention to approximate order-N solutions to (2.71). We begin by fo rma l i zing the notions 
of “a discrete polynomial approximation” and “convergence.” 

Definition 2.1 Let (x„, y„) be a point in E 2 , and let Ax and Ay be positive numbers. 
Let NoaxAv = {(*,y) : |*-*o| < ^ and |y-y 0 | < ^}- Then, alocal discrete polynomial 

approximation to the exact solution of equation (2.71) is a function h N = (J} Z n-, h y N ) 

defined on No^iAy by 

h x s(x,y) = ^ ^ h*_ fe Jfc (Ax, Ay) (x — x 0 ) n * (y — y 0 ) (2.82) 

n =0 fc =0 
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(2.83) 


Tf n 


h'»(z,y) = ^2^2 h n-k,k(^ x ’ A v)( x - x '>) n k (y~y «)*> 

n =0 fc =0 

where h*_ fc Jfc and h^_ kk axe functions of Ax and Ay. 

— * 

Definition 2.2 A local discrete polynomial approximation h N converges to order K 
to the exact solution h of the conservation law V • 7t = 0 as Ax — ► 0, Ay — ► 0 if and 
only if for any e > 0, there exist numbers S\ > 0 and 62 > 0 such that when Ax < and 
Ay < S 2 , 


|k;_ M (A*,Ay) - 

*’-»•*( A*, A,) - if 


< c 


< c 


(2.84) 

(2.85) 


— * — # —# — * 

for n = 0, 1, A, and fc = 0, 1, n, and the remainder R K = h s — h K — ♦ 0 as Ax — ► 0, 
Ay — *■ 0. 

We may now state and prove the two following important theorems. 

Theorem 2.1 Let e be any positive number, and let 0 < K < N. If h w (Ax,Ay) 
converges to order K to the exact solution h of the conservation law V • TC = 0 as 

Ax — ► 0, Ay — * 0, then for all sufficiently small Ax and Ay, 


II £ - ?*(Ax, Ay)|| OO ^ 


Proof: We have 

^ = h hff ft# — hpf (2.2a) 

^ 11 ^ — b#\\oo — ||^ “ || 00 + \\h N — h N || 00 (2.2b) 

^ ||^ “ It 00 ^ ||^ ”* ^n||oo + \\hjc ~ 1 }k || 00 + || Rk “ -?k*||oo ( 2 . 2 c) 

-* -# -♦ -* -♦ -♦ 

where — h N — h K and R K = h N — h K . The first two terms on the right hand side 
of (2.2c) are each less than | for sufficiently small Ax and Ay, and the third term is a 
polynomial whose lowest order term is of degree K 4- 1. Since the coefficients of R K are 
fixed, and R K — ► 0, the third term is also less than | for sufficiently small Ax and Ay, 
and the theorem follows. 

Theorem 2.2 Let e be any positive number, and let 1 < K < N. If h^(Ax, Ay) 
converges to order K to the exact solution h of the conservation law V • Ti = 0 as 

Ax — ► 0, Ay — ► 0, then for all sufficiently small Ax and Ay, 

(n-k) h*_ M (Ax % Ay) + (k + 1) /»*_*_ M+1 (Ax, Ay) < e 
for n = 1 , 2, ...K, k = 0, 1, ...n — 1 . 
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Proof: 

Since h N ( Ax, Ay) converges to order K , for any n = 1, 2, ...,K, and k = 0,1, n — 1, 
we have d n h z 


fcn-k,Jb( Ax ’ A y) 


Ldx"-*dy* (X °’ yo) ] («-*)! k\ 

, . . x r d n h* 0 1 1 

ft n-jfe-i,it+i( Ax ’ A y) “> Lax B - Jfc - 1 3y fc+1 ^ 0 ’ y ° J(n-Jb-l)! (fc + 


so 


that 


1)! (k + 1)! 

(n — k) h*_ fc> *(Ax, Ay) - [ ax „_fc 5y * ( x °’ y °)] ( n _ * - l)! k\ 

{k + 1) h y n _ k _ l k+1 {Ax,Ay) -* [ 5x n-t-i 9 y i+i( Xo ’ y °)] ( n - k - 1)! k\ 
=> (n — k) h^_ k>k ( Ax, Ay) + (k + 1) h*_ k _ l k+1 (Ax, Ay) 

d n h x , x d n h* , /| 1 1 

( x o,yo) + 3 x n— *— I^fc+T^O’H ( n — fc — 1)! Jk!* 


(2.5a) 

(2.56) 

(j8.5c) 

(2.M) 

( 2 . 5 e ) 


Ldx"-*dy* V ~'” i ' 0 ' 1 dx n ~ k ~ 1 dy k+ 

But since h is a solution of the conservation law, 

[ ax lt ay »( x °-!'.) + ax»-^-W+ T(a!, ’ ! '” ) ] = 0 

by Lemma 2.1. Thus, 

(n — A:) fc (Ax, A y) + (* + 1) h* n _ k _ ltk+1 (Ax, Ay) -* 0 


(*tf) 




as Ax — » 0, Ay — » 0, and the theorem follows. 

The meaning of Theorem 2.2 becomes clear when we consider the divergence of h K . 

With h N d = (h x N , h 9 „) defined by (2.82) and (2.83), we have 


d_ 

dx 


d_ 

dy 


N n— 1 

EE<- 

n=l k—0 
N n 

= EE 

n=l t=l 


N n— 1 


n=l fc=0 


so that 


d_ 

dx 


,,*(* - X .)"- 1 * 1 (V - ».)* 

(2.86) 

x - x .)"-‘(»- y .)*-‘ 


- x .)”" 1-1 (y - y.) t 

(2.87) 

S- fc»(x,y) 

oy 

(2.88) 
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= 53 X) [(** - fc ) + (* + *) (* - *«)" * 1 (l/ - Vo)*. 

n=l Jb=0 

According to Theorem 2.2 , a necessary condition for h N to converge to order N is that 

(» - *) K-k,k + (* + 1) />L»- i,.+i - 0 (2-89) 

as Ax, Ay — ► 0, for n = 1, 2, N , k = 0, 1, n — 1. 

The implications of this are especially significant when it comes to numerical calcu- 
lations. In general, the mechanism whereby conditions (2.89) are satisfied depends on the 
particular numerical method being used. Finite-difference methods, for example, satisfy 
each of the conditions (2.89) to a certain order through the difference approximations 
that are used. In this case, condition (2.89) is satisfied to a given order, say order L, for 
n = 1. Then for hi g her values of n, conditions (2.89) are satisfied to an order which is 
less them or equal to L. The higher order constraints expressed by (2.89) do not result in 
independent conditions for a finite-difference scheme. Rather, the higher order constraints 
are automatically satisfied by virtue of the difference equations employed to satisfy (2.89) 
corresponding to n = 1. 

On the other hand, when one solves for the unknown coefficients and h y n _ kk 

directly, as in the present approach, each constraint associated with (2.89) represents an 
independent condition. Thus, to ensure that (2.89) is always satisfied, one should require 

K -k,k h l-k,k to 

(n — k ) + (k + 1) = 0 (2.90) 

for n = 1, 2, ..., N, k = 0, 1, ..., n — 1. That is, hit should be a solution of the conservation 
law, and for the special case of N = 2, constraints (2.68) and (2.69) should always be 
imposed. As a result, conditions (2.89) are not satisfied just to a certain order as in 
finite-difference methods, but rather axe satisfied identically. 

Since the coefficients k and h^_ k k are functions of intermediate variables, each 
constraint associated with (2.90) must be expressed in terms of the intermediate vari- 
ables. For the Navier-Stokes equations, this can be accomplished using (2.25) - (2.27). 
Equations (2.68) - (2.69), corresponding to the conservation of mass, x-momentum, and 


y-momentum, then become 

2uj 0 + v M = 0 (2-91) 

u x>l + 2i? 0> a = 0 (2.92) 

2(2u 0j0 Uj o + + Pa.o) + Ui,i»o,o + Vi,i«o,o + « 1 , 0 » 0,1 + = 0 (2.93) 

2 -b 2u 1>0 « 0>1 -|- Pi,i ~b 2(uo o Vo , 4" "b ^o,i^o,i) 0 (2.94) 

2(u 0i oV 2 ,o + Uo,otta,o + t>i,o Ui.o) + 2v ltl v 0i0 + 2v li0 v 0>l + p M = 0 (2.95) 
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U,,! v 0|0 + Vx.iUo.o + tli.oVo.i + “o,i v i,« + 2(2 V 0 ,« v o ,2 + V* A + Po,i) — 0> (2.96) 

respectively. Note that these constraints ensure the satisfaction of necessary, but not 
sufficient, conditions for the discrete second derivatives to converge to their analytical 
counterparts. 


III. High-Revnolds-Number Flow Past a Finite Flat Plate 

We now apply the discretization outlined above to the thin airfoil boundary layer 
problem. Consider the mesh shown in Figure 3. The airfoil lies on the x-axis between 
x = 0 and x = 1, and has uni form spacing in each direction. Note that the mesh includes 
an upstream and wake region. 

Let Ni and Nj denote the number of solution elements in the x and y directions, 
respectively. Then, representing u, v, and p by way of the local expansions (2.13) - (2.15), 
there are 18 discrete unknowns per cell, for a total of 18 NiNj unknowns altogether. We 
thus require 18 NiNj conditions to have a closed system of equations. 

The flux conservation constraints (2.56), (2.59), and (2.63), together with the differen- 
tial constraints (2.91) - (2.96), immediately provide 9NiNj conditions. These constraints 
ensure that the discrete flux vectors defined by (2.25) — (2.27) satisfy both the integral and 
differential forms of the governing conservation laws. 

To ens ure that mass and momentum are conserved globally, we must require that fluxes 
be conserved across cell interfaces. Using (2.51) and (2.53), the generalized statement of 
flux conservation across vertical interfaces is given by 




Ay 2 ir , Ax 2 Lr Ax fcI ^ h x } _ n 

- -J2 - ^0,2 + “ K,0 ~ 2 hl -° /l °'°J i U ‘ 

For conservation of mass, (3.1) becomes (using (2.55)) 

rA y 2 Ax 2 , Ax ] 

hr-** + + T" 1 ' + -i 


(3.1) 


rAy 2 Ax 2 Ax 1 _ n 

—rz~ Uo,j + — — « 2 ,o _ -T- u i,o + u o,o . — 0* 

12 4 2 J* 


Ax 

2 


(3.2) 


A similar statement is obtained for mass conservation across horizontal interfaces. Using 
(2.58) and (2.62), one derives the analogous expressions for conservation of x and y mo- 
mentum. Thus, one obtains a total of 3 Nj(Ni — 1) + 3 Ni(Nj — 1) — 3 N a interface flux 
conditions, where N a is the number of solution elements between the airfoil leading and 
trailing edges. 

Boundary conditions account for an additional 47V, -f- 3N j + 47V a conditions. For each 
cell adjacent to the airfoil, we require the mass flux through the wall face, and the u velocity 


15 



component at the midpoint of the wall face, to be zero. At the upstream boundary we 
specify the velocity, and at the downstream boundary we specify the pressure. Along the 
free-stream boundary cells, we specify zero y gradient conditions for u and v. 

Finally, as a means of imposing the zero (constant) pressure gradient condition, we 
set 

Pi,i = 0 Pa, o = 0. (3.3a, b) 

These conditions are both justified when the Reynolds number is large, since the pressure 
is nearly constant except near the singular leading and trailing edge regions. A scheme 
which retains p 1A and p 3>0 has also been developed and implemented (Scott et al, 1995). 
Numerical results for the flat plate problem have shown that p 1(1 and p 3|0 are negligible for 
Reynolds numbers greater than about 7500. 

The above conditions ensure the satisfaction of local and global flux conservation, 
boundary conditions, and all other relevant physical requirements. We are now free to 
impose any other physically realistic condition to dose the system. The number of condi- 
tions needed is Ni(Nj — 1 ) — N a . This is precisely the number of horizontal interfaces in 
the mesh (minus those that coincide with the airfoil). Consequently, there is an additional 
degree of freedom in specifying horizontal interface conditions. For the present problem, 
we require the u velocity to be continuous at the midpoint of each horizontal interface. 

By virtue of conditions (3.3), there are 16 unknowns on each solution element. How- 
ever, using the local constraints (2.56), (2.59), (2.63), (2.91), (2.92), and (2.96), six of the 
unknowns may be eliminated in terms of other variables. The total number of unknowns 
that must be explicitly solved for is then 10 NiNj. 

The discrete boundary value problem outlined above is a coupled system of second- 
order polynomial, equations in the unknown coefficients u 0|0 , v 0 ,o, p 0 ,«, etc.. Our experience 
has shown that this system may be solved very efficiently using Newton’s method. 

In what follows, we discuss numerical results for high-Reynolds-number flow (Re = 
100,000) around a finite flat plate. Previous efforts to solve the full Navier-Stokes equations 
for the flat plate problem have usually involved a vorticity-stream function formulation 
(Dennis et al, 1965; McLachlan, 1991). To our knowledge, no solutions to the full Navier- 
Stokes equations have yet been published for Reynolds numbers higher than about 4000. 
McLachlan (1991) has reported some of the numerical difficulties that are encountered 
as the Reynolds number approaches 4000. On the other hand, as we stated above, our 
assumption (3.3) becomes invalid for Reynolds numbers less than about 7500. As a result, 
no direct comparison with vorticity-stream function results can be made at this time. 
However, at the very high Reynolds number of 100,000, we may compare with the analytical 
solutions of Blasius (1908) and Goldstein (1929) to assess the accuracy of our results. 

The Blasius solution may also be used as a guide for constructing a mesh for numerical 
calculations. If we let Re denote the Reynolds number based on the airfoil chord c, then 
the similarity relation (Schlicting, 1979; White, 1974) 
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becomes 

n = VRe v - 

c 

where all quantities except r/ and Re are dimensional. Identifying the free-stream velocity 
U with Uqo, and assuming constant viscosity, equation (3.5) becomes 

" = (3 - 6) 
Corresponding to any 77 , x and y axe related through the equation 



Equation (3.7) allows us to estimate the location of the edge of the boundary layer for any 
fixed We take r/ = 6.0, corresponding to = 0.999999, (White, 1974) to be the edge 
of the C boundary layer. An estimate of the boundary layer thickness 8 X at any § is then 
given by 

8 X = 6.0 

At the trading edge we have the estimate 



When Re = 100,000, one gets 8 t . e . = 0.0268. The free-stream mesh boundaries may then 
be taken to be y = ±0.027. Similarly, at the leading edge, say f = .01, the estimated 
boundary layer thickness is 0.00268 « 0.0027. An estimate for Ay is then 0.0027/N;. e ., 
where Nj e . is the number of cells required to resolve the leading edge boundary layer. If 
Ni. e . = 2, one obtains Ay = .00135. (Recall that the present discretization represents u 
quadratically. See (2.13).) 

Unlike the free-stream boundaries, the upstream and downstream boundaries may be 
located somewhat arbitrarily. For the present study, we take the upstream boundary to 
be x = —0.1, and the downstream boundary to be x = 1.5. 

Newton’s method has proven to be a robust solution technique for the present dis- 
cretization. Our experience has shown that an initial guess of uniform flow is usually 
sufficient to ensure convergence. However, convergence problems have been observed for 
certain mesh spacings and grid aspect ratios. In Figure 4 we show the results of a system- 
atic study which was performed to determine the scheme’s domain of convergence at Re 
= 100,000. For purposes of the study, we limited the number of solution elements in the 
x direction to 400, and in the y direction to 60. This corresponds to minimum x and y 
spacings of ^ = 0 004 and ^ = 0.0009, respectively. The mesh boundaries were held 




17 



fixed at y = ±0.027, x = —0.1, and x = 1.5. The scheme is stable and convergent on the 
boundary and interior of the region shown in Figure 4. The scheme remains stable and 
convergent to the exterior of the dashed portion of the boundary, but becomes unstable 
and/or nonconvergent to the exterior of the solid boundary. (We should point out that 
when we used symmetry to solve for only the upper half flow field, the scheme was un- 
stable for all mesh spacings. Apparently, the instability was caused by the leading edge 
singularity, which invalidated our symmetry boundary condition.) 

Using our earlier estimate Ay = .00135 for the y mesh spacing, along with the fixed 
grid boundaries described above, we obtain a baseline grid by choosing a grid aspect ratio 
r = |f. liking r = .155, we obtain a baseline 185 x 40 grid with Ax = .0087. 

For flat plate flow at high Reynolds numbers, it is well known that the Blasius solution 
provides the correct boundary layer profile over a substantial part of the airfoil. Stewart son 
(1969) and Messiter (1970) have shown that the Blasius solution remains valid to within 
an 0(i2e~t) region near the trailing edge. On the other hand, near the leading edge, the 
Blasius solution is not valid but becomes valid as Re x (Reynolds number based on distance 
from leading edge) becomes large (Jones et al, 1988). 

Figure 5 shows how our Navier-Stokes solution differs from the Blasius solution for a 
series of calculations that were performed on successively finer grids. The coarsest mesh is 
the baseline grid described above. Each succesive grid is refined by adding four cells in the 
y direction and holding r constant at .155. The finest grid is 313 x 68 with Ay = .000794, 
A x = jjjj = .00513, and r = .15485. The figure shows ||u 0 ,o — ub||oo = maXj|u 0 ,o — «b| 
as a function of x for the full length of the airfoil, where ub is the Blasius streamwise 
velocity. The results indicate that our Navier-Stokes solution and the Blasius boundary 
layer solution are in closest agreement at around x = .75, the three-quarter chord point. 

In Figures 6 and 7, using results from the baseline and 313 x 68 grids, we compare u 0>0 
and v 0 ,o with their analytical counterparts at x = .25, .5, and .75. The results obtained 
from the finest mesh clearly show excellent agreement for both u 0>0 and u 0j0 . At x = .75, the 
maximum difference between u 0>0 and u# is .00157. If we take ub to be the exact solution, 
we find the maximum relative error of u 0>0 to be .0032 at x = .75. For the baseline grid, we 
find a maximum absolute error of .00486, and a maximum relative error of .0265. Based 
on these results, we estimate that the fine grid calculation is in error by less them .32%, 
and the baseline grid is in error by less than 2.7%. 

We now examine the rate at which the error (i.e., ||u 0(0 — «b||oo at x = .75) is reduced 
with mesh refinement. Figure 8 shows the effect of refining in the x direction while Ay is 
held constant. Surprisingly, the results show that x grid refinement at constant Ay (i.e., 
increasing the grid aspect ratio) leads to an increase in the error. As Ay decreases, this 
effect becomes stronger. Figure 9 shows a log-log plot of the error versus r. The slopes for 
the four cases are .156, .230, .310, and .390. 

It follows from the above that mesh refinement in the y direction must be done at 
constant r if we wish to determine how the error is reduced with Ay. Figure 10 compares 
the error from 148 x 32 and 296 x 64 grids with r fixed at .15525. The maximum error is 
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.007118 and .001754, respectively, so that the error is reduced by a factor of 4.06. 

Figure 11 shows a log-log plot of the error versus Ay for five grids with r = .15525. 
The fine-grid to coarse- grid slope is 2.02. The four intermediate slopes are, from left to 
right, 1.84, 2.24, 2.35, and 1.71, for an average slope of 2.035. The reduction in error is 
thus second order in Ay. 

The above results imply that the error decreases like r a Ay 2 where, roughly, 0.156 < 
a < 0.390. Mesh refinement at constant Ax confirms that this is the case. Refining a 296 
x 32 grid to a 296 x 64 grid reduces both r and Ay by a factor of two. Taking o = .156, 
one predicts that the error will be reduced by a factor of (2 0-1 ® ) (2 ) = 2 4.457. 

From n um erical calculations we obtain an error reduction factor of 001754 4.511, which 

corresponds to a = 0.174. Similarly, refining a 296 x 56 grid to 296 x 64, one predicts that 
the error will be reduced by a factor ranging from (i) 2 156 = 1.334 to (f) 2 ' 39 = 1*376. 
From n um erical calculations we find that the error is reduced by a factor of 001754 1.381, 

which corresponds to o = 0.417. The above results are summarized in Table 1, along with 
additional results. We find that a assumes an average value of 0.247. The one anomalous 
result, where a = -.395, results from an unusually favorable accuracy on the 313 x 60 grid, 
as is evidenced by the last case shown in the table, where a = .586. 

We now present numerical results for the pressure coefficient, coefficient of friction, 
and wake velocity (from the 313 x 68 grid described above). The trailing edge pressure 
singularity shown in Figure 12 indicates that C p assumes a mi nimum value of -.0149. This 
compares well with the results of Srinivasan et al (1992), who studied flat plate trailing 
edge flow at Re = 100,000 using a reduced Navier-Stokes (RNS) scheme. 

The results shown in Figure 13 indicate that our Navier-Stokes C/ is slightly larger 
than the Blasius C f ( C f (B ) = ^g^) near the trailing edge. This agrees with the findings 
of a number of researchers who have derived second-order corrections to the Blasius drag 
(Kuo, 1953; Imai, 1964; Stewartson, 1969; Messiter, 1970). Near the leading edge, our 
results show a small decrease in Of which continues up to x = .0231. The net result is 
that the increase in C/ near the trailing edge is more than cancelled by the decrease at the 
leading edge, and we obtain a drag coefficient of .004171 which is 99.31% of the Blasius 
value of .004200 = 

In Figure 14 we compare with Goldstein’s (1929) wake solution. The solid lines in 
the figure are our Navier-Stokes solution. They were generated using the full functional 
representation of u on each cell (see 2.13). The diamond-shaped symbols represent the 
numerical solution at the cell center (i.e., u 0>0 ), whereas the squares are Goldstein’s solution. 
(For accuracy reasons we limited the comparison to points where Goldstein’s r\ was less 
than or equal to .75.) The results show that the agreement is not very good until around x 
= 1.15. This is consistent with triple deck theory, which has shown that Goldstein’s wake 
solution is not valid in an 0(Re~i) region near the trailing edge. It is clear, however, 
that our numerical results do not agree nearly as well with the Goldstein wake solution as 
they do with the Blasius flat plate solution. This can be attributed to the fact that the 
Goldstein solution is an asymptotic solution to the boundary layer equations, whereas the 


19 



Blasius solution is an exact solution (admittedly obtained by numerical means). 

We conclude with a discussion of results from nonuniform grids. Since the entire 
discretization presented in Section II remains valid if Ax and Ay are replaced by A x< and 
A j/j, respectively, application of our scheme to variably-spaced grids is straightforward. 

The results shown in Figures 15 and 16 are taken from fine and coarse grids, respec- 
tively. The fine grid is 240 x 40 with 4.1% exponential y stretching and has Ay(wall) 

d = A y w = .0009. The coarse grid is 112 x 20 with 14.7% exponential y stretching and 
A y w = .00135. The solid lines in the figures are obtained using the full functional solution 
for the u velocity, and the diamonds correspond to values of u 0>0 . Comparison is made 
with both the Blasius and Goldstein solutions. The accuracy of the fine grid solution is 
comparable to that of the finest uniform grids considered above, and the accuracy of the 
coarse grid solution is comparable to that of the baseline 185 x 40 mesh. We obtain drag 
coefficients Cp = .004166 from the fine grid and C& = .004124 from the coarse grid. Based 
on our earlier drag results, we estimate that these values are in error by less than .4% and 
1.4%, respectively. Both computations converged in eight Newton iterations, starting from 
uniform flow, to a maximum residual error which is less than 1.2 x 10 -10 . The CPU times 
on a Cray YMP were 580 seconds for the fine grid and 50 seconds for the coarse grid. 
These times would be reduced by about 15% if our code was optimized for computational 
efficiency instead of storage efficiency. 

Finally, we should mention that numerical stability for nonuniform grids depends 
rather strongly on the wall grid aspect ratio, r w = Our experience has shown that 

the domain of convergence in Figure 4 very nearly provides the domain of convergence for 
stretched grids if one replaces r and Ay with r w and Ay w , respectively. 

Concluding Remarks 

The numerical results presented above have clearly shown an accurate resolution of 
the flat plate flow field. We showed in particular that the streamwise velocity is obtained 
to order r a Ay 2 , where a « 0.25. In a future paper we will perform a similar analysis to 
determine the order of accuracy of the discrete derivatives. 

Finally, the effort to generalize the discretization presented in this paper is ongoing. 
We recently presented (Scott et al, 1995) a methodology for retaining all the discrete 
unknowns in the local expansions, and extended the discretization to quadralateral con- 
servation elements of arbitrary shape. Numerical results obtained to date would seem to 
indicate that the present approach offers significant potential to resolve laminar boundary 
layers on relatively coarse grids with only moderate computational effort. 
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Figure 3 Computational grid for a finite flat plate. The airfoil lies 
on the x-axis between x = 0 and x=1. 
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it it , i Grid Aspect Ratio 

||wo,o — ub||oo = maxj|u 0 ,o ~ u b| K 


0.75 r- 


Unstable 



Figure 4 Domain of convergence for a finite flat plate at Re = 100,000. 

The scheme Is convergent on the interior and boundary of the enclosed region, 
and remains convergent to the exterior of the dashed boundary. 



Figure 5 Maximum deviation of Navier-Stokes solution from Blasius solution 
at Re = 100,000 for eight different calculations. Grid aspect ratio = .155. 
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Figure 6a Comparison of predicted u velocity with Blasius solution at 
x = .25, .5, and .75. 185 x 40 Baseline Grid. Re = 100,000. 
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Figure 6 b Comparison of predicted v velocity with Blasius solution at 
x = .25, .5, and .75. 185 x 40 Baseline Grid. Re = 100,000. 
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Figure 7a Comparison of predicted u velocity with Blasius solution at 
x = .25, .5, .75. 313 x 68 Fine Grid. Re = 100,000. 
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Figure 7b Comparison of predicted v velocity with Blasius solution at 
x = .25, .5, .75. 313 x 68 Fine Grid. Re = 100,000. 
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Figure 8 Effect of grid aspect ratio on accuracy for four different 
values of Ay . 
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Figure 9 Effect of grid aspect ratio on accuracy for four 
different values of Ay . The slopes are, from top to 
bottom, .156, .230, .310, and .390, respectively. 
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Figure 10 Comparison of numerical error at x s .75 for 148 x 32 and 296 x 64 
grids with r = .15525. Max error = .007118 and .001754, respectively. 



Figure 1 1 Reduction in error at x = .75 due to mesh refinement 
at constant grid aspect ratio r = .15525. Grid sizes = 148 x 32, 
185 x 40, 221 x 48, 259 x 56, and 296 x 64. 
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296 x 32 

to 

296 x 64 

.007913 

.001754 

= 4.511 

2.174 

0.174 

296x40 

to 

296x64 

.005497 

.001754 

= 3.134 

2.430 

0.430 

296 x 48 

to 

296 x 64 

.003481 

.001754 

= 1.985 

2.383 

0.383 

296 x 56 

to 

296x64 

.002422 

.001754 

= 1.381 

2.416 

0.416 

313 x 36 

to 

313 x 68 

.006341 

.001573 

= 4.031 

2.192 

0.192 

313 x 44 

to 

313 x 68 

.004292 

.001573 

= 2.729 

2.306 

0.306 

313 x 52 

to 

313 x 68 

.002784 

.001573 

= 1.770 

2.128 

0.128 

313 x 60 

to 

313 x 68 

.001923 

.001573 

= 1.223 

1.605 

-0.395 

313 x 52 

to 

313 x 60 

.002784 

.001923 

= 1.448 

2.586 

0.586 


Table 1 Grid refinement at constant Ax. a aV g = >247 
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.25 


Figure 12 Pressure coefficient in the trailing edge region. 
Re = 100,000. 313 x 68 Grid. 
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□ Goldstein 
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Figure 1 4 Comparison of predicted u velocity with Goldstein wake solution 
x s 1 .00256, 1 .04872, 1.1, 1. 15128, 1 .20256, 1 .3, and 1 .5. 313 x 68 Grid. 
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Figure 15a Companion of predicted u velocity with Blasius solution at 
x r .25, .5, .75. 240 x 40 Stretched Grid. Re = 100.000. 


3i 




Figure 1 5b Comparison of predicted u velocity with Goldstein wake solution 
at x = 1.00333, 1.05, 1.10333, 1.15, 1.20333, 1.30333, and 1.49667. 

240 x 40 Stretched Grid. Re = 100,000. 
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